#include "FarFieldRiemann.h"
#include "Attitude.h"
using namespace Eigen;

template<>
InputParameters validParams<FarFieldRiemann>()
{
	InputParameters params = validParams<FarFieldBC>();
	return params;
}

FarFieldRiemann::FarFieldRiemann(const InputParameters & parameters):
	FarFieldBC(parameters)
{
}

void FarFieldRiemann::superSonicInlet()
{
	_cfd_data_neighbor.uh[0] = _rho_inf;
	_cfd_data_neighbor.uh[1] = _rho_inf*_vel_inf(0);
	_cfd_data_neighbor.uh[2] = _rho_inf*_vel_inf(1);
	_cfd_data_neighbor.uh[3] = _rho_inf*_vel_inf(2);
	_cfd_data_neighbor.uh[4] = _p_inf/(_gamma-1) + 0.5*_rho_inf*_vel_inf.norm_sq();

}

void FarFieldRiemann::superSonicOutlet()
{
    _cfd_data_neighbor.uh[0] = _cfd_data.uh[0];
    _cfd_data_neighbor.uh[1] = _cfd_data.uh[1];
    _cfd_data_neighbor.uh[2] = _cfd_data.uh[2];
    _cfd_data_neighbor.uh[3] = _cfd_data.uh[3];
    _cfd_data_neighbor.uh[4] = _cfd_data.uh[4];
}

void FarFieldRiemann::subSonicInlet()
{
	Real cl = _cfd_data.c;
	Real cr = sqrt(_gamma*_p_inf/_rho_inf);
	Real s = _p_inf/pow(_rho_inf, _gamma);
	Real Rp = -_vnr+2.0*cr/(_gamma-1);
	Real Rm = -_vnl-2.0*cl/(_gamma-1);
	Real vnb = -(Rp + Rm)/2.0;
	Real cb = (Rp - Rm)*(_gamma - 1)/4.0;

	Real rho = pow((cb*cb)/(s*_gamma), 1.0/(_gamma - 1));
	RealVectorValue mom = rho*(_vel_inf + (vnb-_vnr)*_normals[_qp]);
	_cfd_data_neighbor.uh[0] = rho;
	_cfd_data_neighbor.uh[1] = mom(0);
	_cfd_data_neighbor.uh[2] = mom(1);
	_cfd_data_neighbor.uh[3] = mom(2);
	_cfd_data_neighbor.uh[4] = cb*cb*rho/_gamma/(_gamma - 1)+ 0.5*mom.norm_sq()/rho;
}

void FarFieldRiemann::subSonicOutlet()
{
	Real cl = _cfd_data.c;
	Real cr = sqrt(_gamma*_p_inf/_rho_inf);
	Real s =_cfd_data.s;
	Real Rp = _vnl + 2*cl/(_gamma - 1);
	Real Rm = _vnr - 2*cr/(_gamma - 1);
	Real vnb = (Rp + Rm) / 2.0;
	Real cb = (Rp - Rm) * (_gamma - 1) / 4.0;

	vnb = (Rp + Rm) / 2.0;
	cb = (Rp - Rm) * (_gamma - 1) / 4.0;

	Real rho = pow((cb * cb) / (s * _gamma), 1.0 / (_gamma - 1));
	RealVectorValue mom = rho*(_cfd_data.vel+(vnb-_vnl)*_normals[_qp]);
	_cfd_data_neighbor.uh[0] = rho;
	_cfd_data_neighbor.uh[1] = mom(0);
	_cfd_data_neighbor.uh[2] = mom(1);
	_cfd_data_neighbor.uh[3] = mom(2);
	_cfd_data_neighbor.uh[4] = cb*cb*rho/_gamma/(_gamma - 1)+ 0.5*mom.norm_sq()/rho;

}

void FarFieldRiemann::boundaryCondition()
{
	_vnl = _cfd_data.vel*_normals[_qp];
	_vnr = _vel_inf*_normals[_qp];

	if(_vnl + _vnr <= 0)
	{
		if (_cfd_data.m > 1)
			superSonicInlet();
		else
			subSonicInlet();
	}
	else
	{
		if (_cfd_data.m > 1)
			superSonicOutlet();
		else
			subSonicOutlet();
	}
}
